In season estimation of EONR
Morteza Abdipourchenarestansofla
2024-02-14
library(sf)## Warning: package 'sf' was built under R version 4.3.2## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUElibrary(rsm)## Warning: package 'rsm' was built under R version 4.3.1library(broom)## Warning: package 'broom' was built under R version 4.3.1library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.2options(scipen = 1000)Read data
base_dir = getwd()
data_dir = paste0(base_dir, "/", "data_maplayers_one_field.csv")
data_wkt <- st_read(data_dir, options = "GEOM_POSSIBLE_NAMES=WKT")## options: GEOM_POSSIBLE_NAMES=WKT
## Reading layer `data_maplayers_one_field' from data source
## `C:\Users\BX2A428\JD\Projects\R\In-season-EONR-estimation-with-smart-farming-data\data_maplayers_one_field.csv'
## using driver `CSV'
## Simple feature collection with 1840 features and 66 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 430277.1 ymin: 5501739 xmax: 430537.1 ymax: 5502069
## CRS: NAmaplayers_sf <- data_wkt[,!(names(data_wkt) == "WKT")]
st_crs(maplayers_sf) <- st_crs("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs")
## Covert datatype to numeric
maplayers_sf[, sapply(maplayers_sf, is.character)] <- lapply(maplayers_sf[, sapply(maplayers_sf, is.character)], as.numeric)## Warning in `[<-.data.frame`(`*tmp*`, , sapply(maplayers_sf, is.character), :
## provided 66 variables to replace 65 variablesmaplayers_df = sf::st_drop_geometry(maplayers_sf)Yield-applied N relationship using Quadratic function
#### Wet-Mass vs applied_rate, N-treatment
ggplot(maplayers_df, aes(x = applied_rate, y=wet_mass_idw)) + geom_point(size=3) +
geom_smooth(method = "lm", formula = y ~ poly(x,2), se=FALSE) +
ggtitle("quadratic regression") + xlab("N-applied, kg/ha") + ylab("harvested grain, kg/ha")
Filter the best identified multispectral indicess during feature selection process using ANOVA and Lasso regression.
candidate_VIs = c(
'cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632',
'cloud.coverage.0_LCI_2020.03.25.10.36.59_7632',
'cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632',
'cloud.coverage.0_CRI550_2020.03.25.10.36.59_7632',
'cloud.coverage.0_EVI_2020.03.25.10.36.59_7632',
'cloud.coverage.0_CARI_2020.03.30.10.37.00_8331',
'cloud.coverage.0_CARI2_2020.03.30.10.37.00_8331',
'cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331',
'cloud.coverage.0_NDRE_2020.03.30.10.37.00_8331',
'cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536',
'cloud.coverage.0_LCI_2020.04.04.10.36.59_1536',
'cloud.coverage.0_CRI550_2020.04.04.10.36.59_1536',
'cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536',
'cloud.coverage.0_NDRE_2020.04.04.10.36.59_1536',
"applied_rate",
"wet_mass_idw")
candidate_VIs_df = subset(maplayers_df, select = candidate_VIs)Modeling with RSM
reduced_model_formula = wet_mass_idw ~ SO(applied_rate,
cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632,
cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632,
cloud.coverage.0_CARI_2020.03.30.10.37.00_8331,
cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331,
cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536,
cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536)
full_model <- rsm(reduced_model_formula, data = candidate_VIs_df)
mdl_summary = summary(full_model)
mdl_summary##
## Call:
## rsm(formula = wet_mass_idw ~ SO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632,
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331,
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536,
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536), data = candidate_VIs_df)
##
## Estimate
## (Intercept) -8665.28562
## applied_rate 277.92723
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 327883.37382
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 -52992.87645
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 48258.42439
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 1813.89118
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -235311.44334
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 5941.49507
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 6156.23351
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 72.03097
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 704.62688
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 11.23289
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -7282.89639
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 21.91666
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 362872.49773
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -583776.08334
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 93856.88086
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -7013590.66180
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -303837.76290
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -7744.49312
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -910.07235
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 51900.43518
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 3275.02536
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 15020.05939
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -3492702.94031
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 12528.83245
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -143638.55208
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 854.88002
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 200675.67714
## applied_rate^2 -2.55804
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2 2697229.06188
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2 -16024.82153
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2 1023909.60102
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2 -671.31560
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2 7523218.75660
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2 721.81252
## Std. Error
## (Intercept) 9348.28152
## applied_rate 156.82393
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 292493.87198
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 13314.31653
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 60808.30099
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 2575.15960
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 310413.34477
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 3552.04074
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 2250.78720
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 93.33038
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 606.04170
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 19.62163
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 2322.32627
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 27.74469
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 172054.07224
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 1102246.45955
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 42395.01989
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 6038832.45311
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 65441.75720
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 43744.44000
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 1500.73728
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 194774.54559
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 2532.93673
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 11098.55669
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 1352880.15990
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 12301.17326
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 47821.72872
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 542.65633
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 63366.54314
## applied_rate^2 0.86505
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2 2691649.96891
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2 3573.73923
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2 153496.48236
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2 287.43781
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2 4018842.95373
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2 496.64474
## t value
## (Intercept) -0.9269
## applied_rate 1.7722
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 1.1210
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 -3.9801
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.7936
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.7044
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -0.7581
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 1.6727
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 2.7351
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.7718
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 1.1627
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.5725
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -3.1360
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.7899
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 2.1091
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -0.5296
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 2.2139
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -1.1614
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -4.6429
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -0.1770
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -0.6064
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.2665
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 1.2930
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 1.3533
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -2.5817
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 1.0185
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -3.0036
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 1.5754
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 3.1669
## applied_rate^2 -2.9571
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2 1.0021
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2 -4.4840
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2 6.6706
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2 -2.3355
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2 1.8720
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2 1.4534
## Pr(>|t|)
## (Intercept) 0.354082
## applied_rate 0.076526
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 0.262440
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.00007160280313
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.427524
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.481287
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.448515
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.094560
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 0.006296
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.440343
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.245117
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.567072
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.001740
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.429666
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.035077
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.596438
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.026963
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.245627
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.00000368410010
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.859497
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.544314
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.789912
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.196185
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.176118
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.009910
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.308574
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.002705
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.115348
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.001566
## applied_rate^2 0.003146
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2 0.316443
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2 0.00000778509314
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2 0.00000000003378
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2 0.019626
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2 0.061370
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2 0.146293
##
## (Intercept)
## applied_rate .
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 ***
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 .
## applied_rate:cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 **
## applied_rate:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632
## applied_rate:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331
## applied_rate:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331
## applied_rate:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 **
## applied_rate:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 *
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 *
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 ***
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI_2020.03.30.10.37.00_8331
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 **
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 **
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536:cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 **
## applied_rate^2 **
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632^2
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632^2 ***
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331^2 ***
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331^2 *
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536^2 .
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536^2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.4066, Adjusted R-squared: 0.3951
## F-statistic: 35.32 on 35 and 1804 DF, p-value: < 0.00000000000000022
##
## Analysis of Variance Table
##
## Response: wet_mass_idw
## Df
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 7
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 21
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 7
## Residuals 1804
## Lack of fit 1804
## Pure error 0
## Sum Sq
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 1024415381
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 530111687
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 158775478
## Residuals 2499956661
## Lack of fit 2499956661
## Pure error 0
## Mean Sq
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 146345054
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 25243414
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 22682211
## Residuals 1385785
## Lack of fit 1385785
## Pure error NaN
## F value
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 105.604
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 18.216
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) 16.368
## Residuals
## Lack of fit NaN
## Pure error
## Pr(>F)
## FO(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) < 0.00000000000000022
## TWI(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) < 0.00000000000000022
## PQ(applied_rate, cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632, cloud.coverage.0_CARI_2020.03.30.10.37.00_8331, cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331, cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536, cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536) < 0.00000000000000022
## Residuals
## Lack of fit NaN
## Pure error
##
## Stationary point of response surface:
## applied_rate
## 68.9526134
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632
## 0.1404002
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632
## 0.4542571
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331
## 0.2210777
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331
## 0.1164771
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536
## 0.1401466
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536
## 1.8876305
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 9605894.552479 1891690.024961 24357.837308 1.302358 -363.129399
## [6] -3736.778354 -289463.272496
##
## $vectors
## [,1] [,2]
## applied_rate -0.0004817841 0.0005765301
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 -0.4412944672 0.6937323329
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 -0.0058733357 0.0703332995
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -0.1644532894 -0.6800575429
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -0.0088768487 0.0060642931
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.8819535609 0.2216808467
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.0160841100 -0.0461617788
## [,3] [,4]
## applied_rate 0.004082907 0.9947852226
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 0.177160547 0.0013894104
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.856472383 0.0250415804
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.342167820 0.0004836887
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 -0.053641044 0.0630152248
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.163038067 0.0007407485
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -0.297514600 0.0761685534
## [,5] [,6]
## applied_rate -0.098126742 -0.02745010
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 0.015727534 0.01007352
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.228642726 0.24255836
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 0.013208515 -0.02810980
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.802205000 -0.59110454
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 0.009986534 -0.01855557
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 0.542253845 0.76796576
## [,7]
## applied_rate 0.001661553
## cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632 -0.540608937
## cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632 0.386955620
## cloud.coverage.0_CARI_2020.03.30.10.37.00_8331 -0.626444871
## cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331 0.009935057
## cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536 -0.382085748
## cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536 -0.139582105Regression diagnostic plots
augmented_data <- augment(full_model)## Warning: The `augment()` method for objects of class `rsm` is not maintained by the broom team, and is only supported through the `lm` tidier method. Please be cautious in interpreting and reporting broom output.
##
## This warning is displayed once per session.# Residuals vs Fitted plot
ggplot(data = augmented_data, aes(.fitted, .resid)) +
geom_point() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
ggtitle("Residuals vs Fitted") +
xlab("Fitted Values") +
ylab("Residuals") + labs(title = "Residuals vs Fitted Plot - Best-model")
# Normal Q-Q plot of residuals
ggplot(data = augmented_data, aes(sample = .std.resid)) +
geom_qq() +
geom_qq_line() +
ggtitle("Normal Q-Q Plot of Residuals") + labs(title = "Normal Q-Q plot of residuals - Best-model")
## Compute model accuracy. RMSE percentage
# Calculate the squared differences between actual and predicted values
squared_errors <- (candidate_VIs_df$wet_mass_idw - full_model$fitted.values)^2
# Calculate the mean of squared errors
mean_squared_error <- mean(squared_errors)
# Calculate the RMSE (take the square root of the mean squared error)
rmse <- sqrt(mean_squared_error)
(rmse/mean(candidate_VIs_df$wet_mass_idw)) * 100## [1] 14.58546Compute N-opt
compute_n_opt = function(trainmodel, first.date.CARI2, first.date.ARI2, second.date.CARI, second.date.CRI550, third.date.CARI2, third.date.CRI700){
mdl_summary = summary(trainmodel)
# collected_terms = mdl_summary$coefficients[c(1:14, 30:36), ]
#######################################################################
### Linear terms plus estimated coefficient for N
# c.0 = mdl_summary$coefficients[c(1, 3:8),]
#######################################################################
first_date_CARI2 = mdl_summary$coefficients[3] * first.date.CARI2
first_date_ARI2 = mdl_summary$coefficients[4] * first.date.ARI2
second_date_CARI = mdl_summary$coefficients[5] * second.date.CARI
second_date_CRI550 = mdl_summary$coefficients[6] * second.date.CRI550
third_date_CARI2 = mdl_summary$coefficients[7] * third.date.CARI2
third_date_CRI700 = mdl_summary$coefficients[8] * third.date.CRI700
c.0 = mdl_summary$coefficients[1] + first_date_CARI2 + first_date_ARI2 + second_date_CARI + second_date_CRI550 + third_date_CARI2 + third_date_CRI700
#######################################################################
### Interaction terms plus estimated coefficient for N
# C.1 = mdl_summary$coefficients[c(2, 9:14)]
#######################################################################
first_date_CARI2_interaction = mdl_summary$coefficients[9] * first.date.CARI2
first_date_ARI2_interaction = mdl_summary$coefficients[10] * first.date.ARI2
second_date_CARI_interaction = mdl_summary$coefficients[11] * second.date.CARI
second_date_CRI550_interaction = mdl_summary$coefficients[12] * second.date.CRI550
third_date_CARI2_interaction = mdl_summary$coefficients[13] * third.date.CARI2
third_date_CRI700_interaction = mdl_summary$coefficients[14] * third.date.CRI700
c.1 = mdl_summary$coefficients[2] + first_date_CARI2_interaction + first_date_ARI2_interaction + second_date_CARI_interaction + second_date_CRI550_interaction + third_date_CARI2_interaction + third_date_CRI700_interaction
c.2 = mdl_summary$coefficients[30]
n_opt = (0.417/0.318) * (1/(2*c.2)) - (c.1/ (2*c.2))
if(class(n_opt)== "data.frame"){
names(n_opt) = "optimum_nitrogen"
names(c.0) = "intercept"
names(c.1) = "interaction_terms"
n_opt$intercept = c.0$intercept
n_opt$interaction = c.1$interaction_terms
}
return(n_opt)
}
#####################################################################################
################### Prediction of EONR ############################################
#####################################################################################
final_vis_predictors = c("cloud.coverage.0_CARI2_2020.03.25.10.36.59_7632",
"cloud.coverage.0_ARI2_2020.03.25.10.36.59_7632",
"cloud.coverage.0_CARI_2020.03.30.10.37.00_8331",
"cloud.coverage.0_CRI550_2020.03.30.10.37.00_8331",
"cloud.coverage.0_CARI2_2020.04.04.10.36.59_1536",
"cloud.coverage.0_CRI700_2020.04.04.10.36.59_1536")
VIs_predictors_df = candidate_VIs_df[, final_vis_predictors]
n_inpute = compute_n_opt(trainmodel = full_model,
first.date.CARI2 = VIs_predictors_df[1], first.date.ARI2 = VIs_predictors_df[2],
second.date.CARI = VIs_predictors_df[3], second.date.CRI550 = VIs_predictors_df[4],
third.date.CARI2 = VIs_predictors_df[5], third.date.CRI700 = VIs_predictors_df[6])
hist(n_inpute$optimum_nitrogen, main = "Economic optimum N rate", xlab = "N-kg/ha")
hist(candidate_VIs_df$applied_rate)
sf_final_data <- subset(maplayers_sf, select = final_vis_predictors)
sf_final_data$optimum_n = n_inpute[, 'optimum_nitrogen']
# Define a color palette with red, yellow, and green
color_palette <- colorRampPalette(c("red", "yellow", "green"))
# Create map view with continuous color scale
mapview::mapview(sf_final_data, zcol = "optimum_n", col.regions = color_palette,
legend = TRUE, legend.style = "classic")Calculate Total N
sf_final_data$optimum_n_total_N = sf_final_data$optimum_n / 0.26
# Define a color palette with red, yellow, and green
color_palette <- colorRampPalette(c("red", "yellow", "green"))
# Create map view with continuous color scale
mapview::mapview(sf_final_data, zcol = "optimum_n_total_N", col.regions = color_palette,
legend = TRUE, legend.style = "classic")Plot two points with different EONR
The location indexes we used for visualization purpose are as follows, 1 and 1800
c_value_0 <- n_inpute$intercept[1]
c_value_1 <-n_inpute$interaction[1]
c_value_2 <- mdl_summary$coefficients[30]
# Single value N
N <- seq(30, 90, length.out = 100)
# Quadratic function
Y <- c_value_0 + c_value_1 * N + c_value_2 * N^2
# Calculate N for the maximum
N_max <- -c_value_1 / (2 * c_value_2)
# Calculate corresponding f(N) for the maximum
Y_max <- c_value_0 + c_value_1 * N_max + c_value_2 * N_max^2
# Plot the curve
plot(N, Y, type = "l", col = "blue", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Function")
# Add a vertical line at the maximum
abline(v = N_max, col = "red", lty = 2)
# Add a legend
legend("bottomright", legend = expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), col = "blue", lty = 1)
# Add text annotations for c_value_0, c_value_1, and c_value_2
text(N[1], Y[1], paste("c[0] =", round(c_value_0, 2)), pos = 4, col = "black")
text(N[50], Y[50], paste("c[1] =", round(c_value_1, 2)), pos = 2, col = "black")
text(N[80], Y[80], paste("c[2] =", round(c_value_2, 2)), pos = 1, col = "black")
# Print the results
cat("N for maximum:", N_max, "\n")## N for maximum: 54.48445cat("Maximum value Y(N):", Y_max, "\n")## Maximum value Y(N): 5804.641c_value_0## [1] -1789.041c_value_1## [1] 278.7468c_value_2## [1] -2.55804c_value_0 <- n_inpute$intercept[1800]
c_value_1 <-n_inpute$interaction[1800]
c_value_2 <- mdl_summary$coefficients[30]
# Single value N
N <- seq(30, 90, length.out = 100)
# Quadratic function
Y <- c_value_0 + c_value_1 * N + c_value_2 * N^2
# Calculate N for the maximum
N_max <- -c_value_1 / (2 * c_value_2)
# Calculate corresponding f(N) for the maximum
Y_max <- c_value_0 + c_value_1 * N_max + c_value_2 * N_max^2
# Plot the curve
plot(N, Y, type = "l", col = "blue", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Function")
# Add a vertical line at the maximum
abline(v = N_max, col = "red", lty = 2)
# Add a legend
legend("bottomright", legend = expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), col = "blue", lty = 1)
# Add text annotations for c_value_0, c_value_1, and c_value_2
text(N[1], Y[1], paste("c[0] =", round(c_value_0, 2)), pos = 4, col = "black")
text(N[50], Y[50], paste("c[1] =", round(c_value_1, 2)), pos = 2, col = "black")
text(N[80], Y[80], paste("c[2] =", round(c_value_2, 2)), pos = 1, col = "black")
# Print the results
cat("N for maximum:", N_max, "\n")## N for maximum: 69.87685cat("Maximum value Y(N):", Y_max, "\n")## Maximum value Y(N): 10389.38c_value_0## [1] -2100.952c_value_1## [1] 357.4955c_value_2## [1] -2.55804Plot a bundle of points.
Picking a few locations in the field and ploting the N-yield response.
Candidate locations are 1, 1800, 100, 200, 56, and 440.
c_0 <- n_inpute$intercept
c_1 <-n_inpute$interaction
c_2 <- mdl_summary$coefficients[30]
c_0_subset = c_0[c(1, 1800, 100, 200, 56, 440)]
c_1_subset = c_1[c(1, 1800, 100, 200, 56, 440)]
c_values_0 <- c_0_subset
c_values_1 <-c_1_subset
c_value_2 <- c_2
# Single value N
N <- seq(0, 100, length.out = 400)
# Create an empty plot with specified ylim
plot(N, rep(NA, length(N)), type = "n", col = "white", xlab = "N kg/ha", ylab = "Y(N) kg/ha", main = "Quadratic Functions", ylim = c(1000, 12000))
# Plot the curves
for (i in 1:length(c_values_0)) {
c_value_0 <- c_values_0[i]
c_value_1 <- c_values_1[i]
# Quadratic function
f <- c_value_0 + c_value_1 * N + c_value_2 * N^2
# Add the curve to the plot
lines(N, f, col = i + 1) # Use different colors for each curve
}
# Add a legend
legend("topright", legend = c(expression(Y(N) == c[0] + c[1] * N + c[2] * N^2), paste("Set ", 2:5)), col = 2:(length(c_values_0) + 1), lty = 1)